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ABSTRACT 

We elucidate the interplay between Newtonian thermal relaxation and numerical dissipation, of several dif- 
ferent origins, in flow simulations of hot extrasolar planet atmospheres. Currently, a large range of Newtonian 
relaxation, or "cooling", times (^10 days to ~1 hour) is used among different models and within a single 
model over the model domain. In this study we demonstrate that a short relaxation time (much less than the 
planetary rotation time) leads to a large amount of unphysical, grid-scale oscillations that contaminate the 
flow field. These oscillations force the use of an excessive amount of artificial viscosity to quench them and 
prevent the simulation from "blowing up". Even if the blow-up is prevented, such simulations can be highly 
inaccurate because they are either severely over-dissipated or under-dissipated, and are best discarded in these 
cases. Other numerical stability and timestep size enhancers (e.g., Robert-Asselin filter or semi-implicit time- 
marching schemes) also produce similar, but less excessive, damping. We present diagnostics procedures to 
choose the "optimal" simulation and discuss implications of our findings for modeling hot extrasolar planet 
atmospheres. 

Subject headings: hydrodynamics — instabilities — methods: numerical — planets and satellites: general — 
turbulence — waves 



1. INTRODUCTION 

There are many studies using a "general circulation model" 
(GCM) to investigate the flow and temperature struc ture of 
close-in extrasolar planet atmospheres (e.g., [Sh owman & 
Guillot||2002j |Cho et al l 1 2003 [ |Cooper & Showman||2005t 
Langton & Lau ghlin|2007| [Cho et al.|2008[ |Dobbs-Dixon &| 
Lin||2008| |S howman et al||2008[ |Menou & Rauscher||20"0"9T 



Thrastarson & Cho|2010[|Rauscher & Menou|2010| l. GCMs 
are advanced numerical models that solve a set of coupled, 
nonlinear partial differential equations for the large-scale mo- 
tions of a shallow fluid on a rotating sphere. In these sophisti- 
cated models, numerous parameters are needed to specify the 
representation of heating and cooling in the atmosphere and 
to stabilize the numerical integration. 

Thus far, not much emphasis has been given to the numeri- 
cal aspects of simulations in the extrasolar planet literature, in 
particular their infl uence on the accuracy of th e model results. 
In an earlier paper, |Thrastarson &"C ho (2010) has investigated 
the sensitivity of initial condition on the extrasolar planet at- 
mosphere flows. In this work, the focus is on another signif- 
icant aspect — the subtle, and not so subtle, interplay between 
numerical and physical parameters. It should be noted that, 
while the discussion is basically numerical, this work is rele- 
vant to both theoretical studies and observations of extrasolar 
planets. 

GCMs usually solv e the hydrostatic primitive equations 
(see, e.g., Salby 1996), which filter sound waves so that only 
two important classes of waves remain — Rossby, or planetary, 
waves (which evolve on slow time scales) and gravity waves 
(which generally evolve on time scales much shorter than the 
Rossby waves). The spatial scales of the two classes of mo- 
tions are generally large and small, respectively. Nonlinear 
advection, which has often been used to define a time scale 
in extrasolar planet work so far, has roughly the same time 
scale as the Rossby waves. Generally, the amplitude of grav- 
ity waves, when averaged over the globe, is very small com- 
pared to that of Rossby waves, and most of the kinetic energy 
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is contained in the large-scale, slow motions. 

A long-standing challenge in GCM theory is finding ways 
to deal with fast waves accurately and efficiently. The fast mo- 
tions not only force small timesteps to be taken (increasing the 
"wall time" of the simulations), they also degrade the fidelity 
with which the equations are solved. Moreover, the very in- 
accuracy often causes the calculation to "blow up" (become 
unstable), preventing any solution at all. With certain types of 
numerical algorithms, such as implicit or semi-implicit time- 
integration schemes, the timestep size restriction can be al- 
leviated. But, artificial viscosity and various filters are still 
required to stabilize the integration in general. 

It is well known that, in conjunction with coarse resolu- 
tion, dissipation and filters can produce results that are se- 
ductively misleading — even to the wary modelers. For exam- 
ple, in the classic H eld-Suarez test for th e dynamical core of 
GCMs for the Earth ( [Held & Suarez|1994| , increasing the res- 
olution generally leads to enhanced equatorward shift of wave 
activity (Wan et al. 2006). The shift becomes more evident in 
the simulations with horizontal resolutions > T85 resolution 
(i.e., 85 sectoral and 85 total modes) so that precise jet posi- 
tions, for example, cannot be ascertained at lower resolutions. 
This is a relatively mild example, but it is telling: the more 
extreme forcing condition for extrasolar planets, it can be ar- 
gued, will lead to larger or more sensitive variations, given 
that the models have been designed and tested for conditions 
appropriate for solar system planets. In this backdrop, even 
inter-comparing different GCMs for extrasolar planet work 
becomes non-trivial. 

In this paper we present and discuss examples of interesting 
behavior when a GCM is stressed to its limits, with what may 
be considered a typical hot, spin-orbit synchronized extrasolar 
planet condition. The implications are broad in the sense that 
the lessons are not just limited to studies using GCMs, but 
also other types of global circulation models. The issues are 
present in all of them. 

The basic plan of the paper is as follows. In Section [2] we 
describe the GCM model we use and its setup for the simula- 
tions described in this work. In Section [3~T1 we focus on the 
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interaction between artificial viscosity and the thermal relax- 
ation time, which is an important parameter in the represen- 
tation of t herm al forcing commonly used in current studies. 
In Section [372] we examine sensitivity of the simulations to 
the Robert- Asselin filter, which is used to stabilize the time- 
marching scheme. We conclude in Section [4] summarizing 
this work and discussing its implications for extrasolar planet 
circulation modeling. 

2. METHOD 

2. 1 . Governing Equations 



In this work , we solve the same equations as in Thrastarson 
|& Cho| ( |2010| l. Here we briefly summarize the relevant as- 
pects for the reader. The horizontal momentum equations are 
solved in the vorticity-divergence form: 



dt 
OS 

dt 



k Vxn+D, 



= V-n+V 2 (E + $) + T> s , 



(1) 



(2) 



where £ = £k=Vxvis the vorticity, 5 = V • v is the diver- 
gence, v is the horizontal velocity, k is the vertical unit vector, 
E = (t-t)/2, <& is the geopotential, and 



n = -((+/)kxv-i) 



d\ RT 

di] p 



where / is the Coriolis parameter, p is pressure, T is tempera- 
ture and R is the specific gas constant. The vertical coordinate 
is a generalized pressure coordinate, rj = rj(p,p s ), with p s the 
bottom surface pressure, and r) = Dr)/Dt with 

D d _ d 
— = — + v V + n— 
Dt dt 'drj 

the material derivative. Hydrostatic equilibrium is assumed: 

<9$ _ RT dp 

dr/ p drj ' 



(3) 



and the ideal gas law, p = pRT, where p is density, is taken 
as the equation of state. The mass continuity equation is in- 
tegrated from the bottom (77 = 1) to the top surface, rj t , using 
the boundary conditions ?) = at both the top and the bottom, 
which yields an evolution equation for p s : 



djh 
dt 



dp 
dr] 



drj 



(4) 



Integration of the continuity equation from ?y t to rj yields a 
diagnostic equation for r): 



.dp dp r> (dp 



drj. 



drj dt J \dr] 
The diagnostic equation for u = Dp/Dt is then: 

dp 



■Vp- 



Finally, the energy equation is 
Dr u 
Dt pc p 



Qnet 



drj 



dr]. 



(5) 



(6) 



(7) 



where c p is the specific heat at constant pressure. In the fi- 
nal formulation of the equations, terms involving v are rep- 
resented in terms of the transformed velocity \coscj), where 



(f> is latitude, in order to avoid discontinuities at the poles. 
Also, the equations are formulated using \n(p s ) instead of p s 
to avoid aliasing problems. The T> terms in the vorticity, di- 
vergence and energy equations represent horizontal diffusion, 
discussed in the next subsection. 

2.2. Numerical Algorithm 

To solve the equations described in the preceding subsec- 
tion, we use the Com munity Atmospher e Model (CAM 3.0), 
describ ed in |Collins et al.| (|2004 ) and Thrastarson & Cho| 
( |2010[ ). CAM is a well-tested, highly-accurate hydro dynam- 
ics model employing the pseudospectral algorithm (Orszag 
[T970l|Eliasenetal.|1970| l. 

For problems not involving sharp discontinuities (e.g., 
shocks or, in atmospheric dynamics problems, fronts) and ir- 
regular geometry, the pseudospectral method is superior to 
the st andard grid and particle methods (e.g., |Canuto et al.| 
1988). To equal the accuracy of the pseudospectral method 
for a problem solved with the computational domain decom- 
posed into N grid points, one would need a N th -order finite 
difference or finite element method with an error of 0{NAx), 
where Ax is the grid spacing and O(-) is the asymptotic order 
(e.g., |Nayfeh|1973) . This is because as N increases, the pseu- 
dospectral method benefits in two ways. First, Ax becomes 
smaller, which would cause the error to rapidly decrease even 
if the order of the method were fixed. However, unlike finite 
difference and finite element methods, the order is not fixed: 
when N is doubled to 2N, the error becomes 0[(Ax) 2N ] in 
terms of the new, smaller Ax. Since Ax is 0(1 /N), the error 
for the pseudospectral method is 0[(l /N) N ]. 

Significantly, the error decreases faster than any finite 
power of N since the power in the error formula is always 
increasing as well, giving an "infinite order" or "exponential" 
convergence. This advantage is particularly important when 
many decimal places of accuracy or high resolution is needed. 
Note that in the vertical direction CAM uses a finite differenc- 
ing scheme, as in most GCMs. 

For the spherical geometry, the horizontal representation of 
an arbitrary scalar quantity £ consists of a truncated series of 
spherical harmonics, 

M N(m) 

ax,p) = £]rec(^ im \ 

N(m) n=\m\ 

where M is the highest Fourier (sectoral) wavenumber in- 
cluded in the east-west representation; N(m), which can be 
a function of the Fourier wavenumber m, is the highest degree 
of the associated Legendre functions P"; A is the longitude; 
and, p = sin </>. The spherical harmonic functions, 



Y;;\X,p) = P' n "(ri)e 



(8) 



used in the spectral expansion are the eigenfunctions of the 
Laplacian operator in spherical coordinates: 



n(n + \) 



(9) 



where 



d_ 
Oft 



d_ 
d/i 



d^_ 

p? dX 2 



and R p is the planetary radius. The set, {T„ m }, constitutes 
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a complete and orthogonal expansion basis (Byron & Fuller 
[19921 . 

In the Navier-Stokes equations, the diffus ion terms appea r 
as the Laplacian of the dynamical variables (Batchelor 1967 ). 
In our case, the diffusion i s generalized to the fo llowing "hy- 
perdissipation" form (e.g., Cho & Porvani|1996" l: 



T>x = v 2p [(-l) p+1 V 2p +C 



(10) 



where \ = {C> T} and C = (2/Rfy p is a correction term added 
to the vorticity and divergence equations to prevent damping 
of uniform rotations for angular momentum conservation. In 
the above form, the p = 2 case is sometimes referred to as 
superdissipation. Hyperdiffusion is added in each layer to 
prevent accumulation of power on the small, poorly-resolved 
s cales and to stabilize the integration. 

|Cho & Polvani| ( |1996[ ) describes the effects of various hy- 
perviscositie s (i .e ., different values of p). As discussed in that 
work, a rational procedure for estimating roughly the value of 
^2p can be obtained in the following way. To damp oscilla- 
tions at the smallest resolved scale (set by the truncation wave 
number, n t ), by an e-folding factor in time r^, one requires 
that 



vip = 0< 



1 

Td 



R, 



n,{n t + Y) 



(11) 



Thereafter, the optimal value of z^p is obtained by comput- 
ing the kinetic energy spectrum (see Section|3]l. Note that the 
precise value is problem specific, and the procedure just de- 
scribed should be performed for each problem — as has been 
done in this work. 

In numerical solutions of time-dependent equations, there 
are two main ways of marching in time. Explicit methods 
give the solution at the next time level in terms of an explicit 
expression which can be evaluated by using the solution at 
the previous timestep. Implicit methods, on the other hand, 
require solving a boundary value problem at each timestep. 
Explicit time differencing is a more straightforward numeri- 
cal approximation to the equations. In our model, the time- 
marching is effected using a semi-implicit scheme, a mixture 
of the two methods commonly used in GCMs. In this scheme, 
the equations are split into nonlinear and linear terms, sym- 
bolically written: 



Ik 



= Afm + £(*), 



(12) 



where Af(^f) and C(^>) denote the nonlinear and linear terms, 
respectively, and ^ is the state of a variable in \ = {C> 8, T}- 

For the nonlinear terms, an explicit leapfrog scheme is used. 
This is a second-order, three-time-level scheme. Because a 
second-order method is applied to solve a differential equa- 
tion which is first-order in time, an unphysical computational 
mode is admitted, in addition to the physical one. In simu- 
lations containing nonlinear waves, the computational mode 
can amplify ov er time, generating a time splitting instability 
( |Durran|[T999l l. Robert (1966) and Asselin (1972) designed 
a filter to suppress the computational mode — hence the time 
splitting instability. This filter is applied in the GCM used in 
the present work. It is applied at each timestep so that 



(13) 



where ^ n = 'J'(nAf), an overbar refers to the filtered state, 
and e specifies the strength of the filter. The filter results in 



strong damping of the amplitude of the spurious computa- 
tional mode. However, it also introduces a second-order error 
in the amplitude of the physical mo de with high values of e, 
as we discuss further in Section [3~2l 

Some parts of the equations can be solved implicitly with 
advantage. In particular, the linear parts that produce fast 
gravity waves are treated implicitly in many GCMs, includ- 
ing the one used in this work. This treatment allows a larger 
timestep to be used, as mentioned in SectionfT] H owever , it is 
also at the cost of degraded accuracy (e.g., Durran 1999). 

As can be seen, time-integration of the primitive equations 
is not a straightforward matter, even with a relatively sim- 
ple method like the leapfrog scheme. The theoretical analysis 
of the scheme is equally complex. The stability of the com- 
bined, semi-imp licit leapfrog scheme has been examined by 
Simmons et al.| ( |1978| l, particularly with respect to the basic 
state temperature profile. They find the isothermal basic state 
distribution to be more stable than a spatially-varying distribu- 
tion, with the stability generally increasing with higher basic 
temperature. In the present work an isothermal basic state of 
1400 K is used. 

2.3. Calculation Setup 

In addition to tuneable parameters associated with the nu- 
merical scheme, such as the ones mentioned in the preceding 
subsection, the representations of physical processes also re- 
quire specification of parameters. Many of these are as yet 
poorly constrained by observations or unobtainable fr om first 
principles (see, e .g., d iscussions in Cho et aL] ( 2008| l, |Show- 
man et al. (2008 1, and |Cho| ( |2008| l)[ One example is thermal 
ind c 



forcing (i.e., heating and cooling) due to the irradiation from 
the host star and radiative processes in the planetary atmo- 
sphere, which is represented in an idealized way currently in 
all extrasolar planet atmosphere simulations. Many crudely 
repres ent the forcing by Newtonian relaxation, as in this work 
(e.g.,|Cooper & Showman 2005 ; Langt on~& Laughiin"]|2007| 
[ Showman et al.|2008[ |Menou & Rauscher|20091 |Rauscher &\ 
|Menou|2010||Tnrastarson & Cho|201^ ~ 

In this representation, the net heating term in equation |7]i 
is represented by 



<?net 
Cn 



Tth 



(T-T e ), 



(14) 



where T e = T e (\,(j),r),t) is the "equilibrium" temperature dis- 
tribution and Tfl, is the thermal relaxation (drag or "cooling") 
time. The appropriate values to use for this relaxation time (as 
well as the equilibrium temperature distribution) are poorly 
known and a large range of values has been used in the extra- 
solar planet literature. In several studies, very short relaxation 
times — even less th an an hour — and large T e gradients have 
been used (e.g., Showman et al. |2008 
[20T0l[Thrastarson & Cho|2010| >. T 



Rauscher & Menou 
his represents a rather "vi- 



olent" forcing on the flow, depending on the initial condition. 

In this work, both and T e are prescribed an d barotropic 
(i.e., d/drj = ) and steady (i.e., d/dt = 0). As in |Thrastarson| 
|& Cho| ( |20T0l >, 



71 = T„, + ATlcoscicos A. 



(15) 



with T m = (7d + 7jv)/2 and AT e = (Td-Tn)/2, where 7b and 
Tn are the maximum and minimum temperatures at the day 
and night sides, respectively. All the simulations described 
in this paper have T D = 1900 K, T N = 900 K. Other physical 
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TABLE 1 
Physical Parameters 



Planetary rotation rate 


Q 


2.1 xlO" 5 


s" 1 


Planetary radius 


R P 


10 s 


m 


Gravity 


g 


10 


m s~ 2 


Specific heat at constant pressure 


c p 


1.23X10 4 


J kg" 1 Kr 1 


Specific gas constant 


R 


3.5xl0 3 


J kg" 1 K" 1 


Mean equilibrium temperature 


T m 


1400 


K 


Equilibrium substellar temperature 


T D 


1900 


K 


Equilibrium antistellar temperature 


r N 


900 


K 


Initial temperature 


T 


1400 


K 



parameters chosen are based on the close-in extrasolar planet, 
HD209458b (see Table [I). 

The spectral resolutions in the horizontal direction for the 
runs described in the paper are T85 and T21. The number 
refers to the maximum total wavenumber, n, = max{N{m)}, 
at which expansion ([SJl is truncated (e.g., T85 => n t = 85); 
"T" means the truncation is such that M = N in equation 
a "triangular truncation" in wavenumber space. A T85 spec- 
tral resolution corresponds roughly to 800x400 grid points in 
physical space of grid-based methods (roughly 200x 100 for 
T21 resolution). The vertical direction is resolved by 26 cou- 
pled layers, with the top level of the model located at 3 mbar. 
The pressure at the bottom r\ boundary is initially 1 bar, but 
the value of the pressure changes in time. The entire do- 
main is initialized with an isothermal temperature distribu- 
tion, To = T m = 1400 K. The flow field is initialized with a 
small, random perturbation; specifically, values of the v com- 
ponents are drawn from a Gaussian random distribution cen- 
tered on zero with a standard deviation of 0.05 m s" 1 . The 
sensitivity to initial flow is described in detail in Thra starsori] 
|& Cho| ( |20T0l l. 

3. RESULTS 
3.1. Spatial Dissipation 

Table |2]lists all the runs discussed in this subsection. Sim- 
ulations are performed with the setup described above, but 
with varying strength of artificial viscosity ( v and p) and the 
forcing timescale {t±). Figure [T] presents the relative vortic- 
ity field near the p sa 85 mb level; this is approximately a 
quarter of the way down from the top of the computational 
domain. The field at t = 80 t p is shown in cylindrical equidis- 
tant projection, centered at the equator, for ten simulations in 
which the setup is identical except for the values of r t h and v 
(p = 2, superdissipation); here, t p = 2ir/Q, is the planetary ro- 
tation period. Positive vorticity (red color) signifies local ro- 
tation in the same direction as the planetary rotation (counter- 
clockwise in the northern hemisphere), and opposite for the 
negative vorticity. The panels on the left column all have the 
same short value of t± = 0.1 t p , while the panels on the right 
column all have t± = 3 t p for five different values of v. In 
all the runs shown, the global kinetic energy time series have 
reached stationary ("equilibrated") state and do not change 
qualitatively for approximately 300 t p . 

For a given value of r t h, simulations with different z/s gen- 
erally share some common features over a range of z/s. But, 
there are clear differences in the character of the flow and 
temperature fields. The differences, which are both qualita- 
tive and quantitative, arise from the strength of dissipation. 



TABLE 2 
List of Runs Discussed 



Run is 2p [m 4 s '] p m/r,, 



Nla 


lxlO 24 


2 


.1 


Nib 


lxlO 24 


2 


3 


N2a 


lxlO 23 


2 


.1 


N2b 


1 x 10 23 


2 


3 


N3a 


1 x 10 22 


2 


.1 


N3b 


lxlO 22 


2 


3 


N4a 


lxlO 21 


2 


.1 


N4b 


lxlO 21 


2 


3 


N5a 


1 x 10 20 


2 


.1 


N5b 


1 x 10 20 


2 


3 


N6a 


t6xl0 12 


1 


.1 



NOTE. — v is the hyperviscosity 
coefficient and p the order index of 
the hyperviscosity. t± is the ther- 
mal relaxation timescale and r ; , is 
one planet rotation. All the simula- 
tions are run at T85 resolution with 
a timestep At of 60 s and a Robert- 
Asselin filter coefficient e of 0.06. 
tThe units for this u are [m 2 s -1 ]. 

Moreover, v can affect the temporal behavior as well. For ex- 
ample, temporal variability can be muted with larger v. Not 
surprisingly, in the strongest dissipation cases [panels (a) and 
(b)] variability in time is essentially completely quenched and 
the flow structures are quite smooth in appearance. These are 
examples of runs which are severely over-dissipated. 

At the other extreme, runs can also be severely under- 
dissipated. This is shown in panels (i) and (j) in Figure [T] 
Note that the common v value in these runs is four orders of 
magnitude smaller than that for the runs of panels (a) and (b). 
A quick visual check of panels (i) and (j) immediately shows 
the physical fields dominated by small-scale oscillations: this 
is numerical noise. Here, by "small" we mean scales near the 
grid-scale, / = O(jSx). Typically , run s like these blow up — or 
at least they should (see Section [3~2| ). Simulations often blow 
up long before the small-scales contain any significant amount 
of energy compared to the large-scales. As we discuss more 
later, this is because the calculation correctly becomes unsta- 
ble. But, sometimes misbehaving simulations can be surpris- 
ingly resilient and not crash. This is usually a signal that bad 
numerics is at play. 

As expected, increasing v leads to decreasing small-scale 
oscillations and to increasingly smoother fields. However, 
significantly, we note that for a given value of v for the two 
Tth's (cf., panels of the same row in Figure [TJ shorter t± in 
a run admits much more pronounced grid-scale oscillations. 
For example, with v\ = 10 22 m 4 s" 1 [panels (e) and (f)], the 
viscosity is clearly insufficient to suppress small-scale oscil- 
lations in the case of = 0.1 r p , while no small-scale oscil- 
lations are present in the calculation with longer r t h. More 
importantly, a v value which appears to be acceptable for 
the shorter relaxation time [e.g., panel (a)] is clearly over- 
dissipative for the run with the longer r t h [e.g., panel (b)]: 
here, the calculation in (b) should be compared with that in 
(f), which is clearly a much less dissipated run than that in 
(b). Hence, running a simulation at a single t± — even if v 
were varied — would not produce trustworthy results since the 
parameter space is at least two-dimensional. 

The implication of this is serious. In many current sim- 
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FIG. 1. — Vorticity field at t = 80t p (planet rotations), near the p « 85 mb level, for two sets of five simulations (left column and right column) that are set 
up identically, except for the viscosity coefficient and the thermal relaxation time. The global kinetic energy time series for all the runs have reached stationary 
("equilibrated") states and do not qualitatively change for ~ 300 r p . The superdiffusion coefficient is z>4 = {10 24 , 10 23 , 10 22 , 10 21 , 10 20 } m 4 s'\ decreasing from 
top to bottom in each column. The panels in the left column have a relaxation time of T t t, = 0.1 t p while the panels on the right have rn = 3 r p . Red (blue) color 
represents positive (negative) values of vorticity, with units s . 
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ulations of hot planet atmospheric flows, a range of r t h's is 
specified, spread over the model atmosphere domain, which 
always contains a region with a short rn,. This forces those 
model atmosphere calculations to be excessively noisy and 
excessively dissipated, in different atmospheric regions of the 
computational domain. Once noise appears in the calculation 
somewhere in the domain, the entire domain becomes quickly 
contaminated. Note that an inherently smooth field — such as 
temperature compared to vorticity, for example — would not 
reveal the noise as well, since it is essentially two integra- 
tions (a smoothing operation) of the vorticity field. In other 
words, temperature possesses a steep (narrow) spectrum like 
the stream function, as opposed to a shallow (broad) spec- 
trum like the vorticity. Similarly, other averaging (integrating) 
procedures, such as taking zonal (eastward) and/or temporal 
means, would obscure, possibly mislead, the analysis of the 
simulation if unaveraged "higher-order" fields like vorticity 
are not considered concomitantly. 

We wish to emphasize here that, in contrast to what might 
be the customary view, numerical noise and blow-ups are use- 
ful. Simulations with severe forcing should be allowed to 
crash — or at least halted, when near grid-scale oscillations are 
visible in the flow field. Any phenomena observed thereafter 
would be seriously compromi sed in a ccuracy, and quite pos- 
sibly entirely artifactual (Boyd 2000). In numerical work, it 
is easy to get lured into believing a calculation by not heeding 
important telltale signs. 

There is a rational way to diagnose the onset of the small- 
scale error sources — as well as the excessive dissipation — in a 
simulation. This is illustrated in Figure [2] which contains the 
kinetic energy spectra of the fields presented in Figure [TJ To 
the best of our knowledge, this is the first time kinetic energy 
spectra are shown in extrasolar planet atmosphere flow sim- 
ulations. They provide an important diagnostic, whe n used 
in conjunction wit h instantaneous fields (see, e.g., |Cho &| 
|Polvani| ( |1996[ l and Koshyk et al.| ( |1999[ l for a discussion of 
kinetic energy spectra and horizontal diffusion), and can be 
used to choose an appropriate v value. 

The left set of spectra in Figure [^corresponds to runs with 
the shorter r t h = 0.1 in the left column of Figure [TJ and 
the right set of spectra in Figure [2] corresponds to runs with 
the longer r t h = 3 r p in the right column of Figure [TJ Vi- 
sual inspection of the vorticity fields along the left column of 
Figure [TJ suggests the runs in panels (a) and (c) are not much 
affected by the small-scale oscillations [if at all in the run of 
panel (a)]. This can be quantified by confirming that the corre- 
sponding spectra in Figure[2](left panel) are the blue and black 
lines (runs Nla and N3a, respectively). In fact, the blue line 
clearly reveals a case of over-damping, in which all scales are 
less energetic than the corresponding scales in the other runs. 

In contrast, note the appearance of near-grid-scale waves in 
physical space, for the run in panel (i) in Figure [TJ indicated 
by a tendency for the spectrum (red line in left panel of Fig- 
ure [2]l to peel off and curl up near — and considerably to the 
left (larger scale) of — the aliasing limit; this is 

2ttR„ 



3Ax ' 



which is 85 in our case, s ince Ax is chosen to be "alias- 
free" up to n t ([Orszag 197TJ. Clearly, our de-aliasing pro- 
cedure, of inverse transforming onto a physical grid that is 
3n t + 1 around the longitude, is not successful in runs N5a and 
N4a, as well as in run N3a (spectrum not shown). This is be- 
cause increasingly greater resolution is needed as the calcula- 



tion proceeds, as discussed below. In turbulence simulations, 
this peeling off behavior is known as an "energy pile-up" or 
"spectral blocking" (because direct energy cascade to high 
wavenumbers in three-dimensional turbulence is blocked). It 
is not limited to spectral methods. It is universal to all meth- 
ods which discretize space. 

Spectral blocking can cause numerical instability in the 
time integration of any nonlinear equations. The instability 
arises due to the quadratically nonlinear term in the solved 
equations. For example, a typical quadratically nonlinear 
term (in one-dimensional Cartesian geometry for simplicity) 



gives: 




2A~ 



k=-2K 



ikx 



Here, ip(x,t) is an arbitrary one-dimensional scalar function, 
which is Fourier expanded; bk are given by a sum over the 
products of the a^, K is the truncation wavenumber, corre- 
sponding to n, in equation ( [TT| . Note that the nonlinear in- 
teraction generates high wavenumbers, k > K, which will be 
aliased into wavenumbers on the range k g [— K,K], This 
induces an unphysical inverse cascade of energy from high 
wavenumbers to low wavenumbers. 

It is important to realize that the above cascade injects ar- 
tificial energy into all scales. The injection is simply more 
noticeable in the small scales since not much energy is con- 
tained there in the absence of blocking. Oscillations of size 
I = O(Ax) are a precursor to breakdown of computational fi- 
delity. These oscillations are insidious because they require 
higher and higher resolution in the calculation over time. 
Without the increasing resolution, they deteriorate the accu- 
racy of the simulation o n all scales as the calcu lation pro- 
ceeds, as pointed out in Thrastarso n & Cho| (f2010). Although 
some blocking is almost inevitable in a long time integration 
of a nonlinear system (unless the dissipation is unrealistically 
large), it can be monitored and controlled — albeit better in 
some methods than in others. 

The left and right panels of Figure [2] reveal not only how 
the appropriate dissipation can be chosen, but also the crucial 
interplay between the small-scale noise and r t h — hence, un- 
derscoring the importance of using both the spectra and the 
physical field in analyzing a calculation. Consider, for exam- 
ple, the "optimal" calculation (i.e., least affected by too much 
or too little dissipation) for the short t± runs. The calculation 
with i>4 = 10 23 m 4 s -1 (black line in the left panel) is devoid of 
non-physical build up of energy at the smallest scales while 
still retaining the same amount of energy in the large scales 
as in the calculations with smaller v. On the other hand, the 
calculation with longer r t h but same v (black line in the right 
panel) is clearly over-dissipated, containing less energy com- 
pared to the other calculations on essentially all the scales. 
Hence, if the v value were "tuned" with the calculations with 
shorter rth (only), then a calculation with a different t± (say 
a longer one, as in this example) would be over-damped. In 
other words, a correct v value cannot be obtained independent 
of r t h. Actually, v = v (r t h, £,•••)> where "• • • " includes T e , R p . 
semi-implicitness, etc. 

The above behavior is generic. Simulations performed with 
a greater range of t± (down to 0.01 t p ) and v and p, exhibit the 
same basic behavior; and, it is present throughout the model 
domain; grid-scale oscillations can appear in the duration of 
a calculation anywhere in the domain. These oscillations can 
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FIG. 2. — Kinetic energy spectra for the fields shown in Figure[T]for simulations that are set up identically, except for the artificial viscosity and the thermal 
relaxation time. The runs shown in the left panel have a relaxation time of = 0.1 Tp while the runs on the right panel have r,), = 3 r p . The different colored lines 
are for different values of v, as indicated in the legend. Note the different scales on the two panels — much more kinetic energy is contained in the flow when the 
relaxation time is short. The spectra reveal both under-dissipated (e.g. red line, left panel) and over-dissipated (e.g. blue line, left panel) flow fields. A v-value 
that seems to give a reasonable spectrum for the short (e.g. black line, left panel) results in over-dissipation for the longer r t i, (black line, right panel). 




1 10 100 

Wavenumber, n 



FIG. 3. — Kinetic energy spectra for simulations that are set up identically, 
except for the artificial viscosity. The different lines refer to different values 
of v, as indicated in the legend. The blue and black lines are the same as 
in the left panel of Figure [2] for which the viscosity is of biharmonic form 
(V 4 with p = 2). But, the red line is for a simulation where the order of the 
viscosity is lower (p = 1), the normal Newtonian viscosity. 



be controlled to some degree in mild cases, as outlined above. 
However, grid-scale oscillations are dominant near the top of 
the domain for all values of v considered. In this situation, 
it is common in GCM studies to include a "sponge layer", 
where dissipation is artificially enhanced in the topmost lay- 
ers. While this can damp unphysical oscillations, it can also 
have spurious effects. The effects of "sponges" as well as 
other boundary conditions will be described elsewhere. 

Figure[3]shows how the spectrum is affected when the form 
of the artificial viscosity is of lower order. The blue and black 
lines (runs Nla and N2a, respectively) are the same as in the 
left panel of figure Figure [2] They can be compared to the 
red line (run N6a), which shows the spectrum from a sim- 
ulation that is identical to the other two runs in the figure, 
except for the value of v and the order of the viscosity oper- 
ator (here p = 1). In this case, the energy in the small scales 
(high wavenumbers) is dissipated much more strongly. More 
importantly, essentially all wavenumbers are affected by the 
lower order viscosity; and, as discussed in |Cho & Polvani| 
( |1996[ ), the slope of the spectrum becomes steeper — even at 
wavenumbers well below the truncation scale. 

3.2. Temporal Dissipation 

If the solved equations support several types of waves, as 
with the primitive equations, the maximum stable timestep is 
limited by the Courant number, 

*= (*\ 



where c max is the maximum horizontal wind speed associated 
with the fastest propagating wave. Some fast waves are of lit- 
tle physical significance, but they enslave Af to be small. Im- 



8 



THRASTARSON & CHO 



TABLE 3 
Summary of Runs for Time Filter 
Sensitivity 

Notes 



Run 


c 


r th/ 


Ela 


0.001 


.1 


Elb 


0.001 


3 


E2a 


0.002 


.1 


E2b 


0.002 


3 


E3a 


0.006 


.1 


E3b 


0.006 


3 


E4a 


0.01 


.1 


E4b 


0.01 


3 


E5a 


0.06 


.1 


E5b 


0.06 


3 


E6a 


0.1 


.1 


E6b 


0.1 


3 



blow-up (f = 1 T p ) 
blow-up (f = 71 T p ) 

bloW-Up (f = 1 Tp) 

blow-up (f = 8 Tp) 



NOTE. — T t (,/r p is the thermal relaxation 
time in units of planetary rotations, and e 
is the Robert-Asselin filter coefficient. All 
the runs are at T21 resolution and have v$ — 
10 22 m 4 s -1 . The timestep is Af = 240 s. 



plicit schemes do permit a larger timestep size to be used than 
in explicit schemes, often making the former more computa- 
tionally efficient. However, for nonlinear equations, implicit 
schemes have a high cost per timestep because a nonlinear 
boundary value problem must be solved at each timestep. 

As noted, a semi-implicit algorithm is commonly used in 
GCMs. In general, the implicit and explicit parts in the al- 
gorithm may be of same or different order. Treating some 
terms explicitly while others implicitly may appear strange, 
but there are some major advantages. First, because the non- 
linear terms are treated explicitly, it is only necessary to solve 
a linear boundary value problem at each timestep. Second, the 
hyperdissipation terms, which involve even number of deriva- 
tives, impose a much stiffer timestep requirement than the ad- 
vective terms; for example, Af is (9(1 /N 4 ) and 0(1 /N 2 ), re- 
spectively, for the Newtonian viscosity (p = 1). Hence, the 
semi-implicit algorithm stabilizes the most unwieldy terms. 
Third, in general circulation and other fluid dynamics prob- 
lems, advection is crucial; therefore, it is important to use a 
high order time-marching scheme with a short timestep to ac- 
curately compute phenonmena or structures such as fronto- 
genesis, advection of storm systems, and turbulent cascades. 
There is little advantage in treating the nonlinear terms im- 
plicitly because a timestep longer than the explicit advective 
stability limit would be too inaccurate. 

Note that, although it is possible to treat the time coordinate 
spectrally, it is generally more efficient to apply spectral meth- 
ods to the spatial coordinates only because time marching is 
usually much cheaper than computing the solution simultane- 
ously over all space-time. In general, much less concern is 
given to the temporal accuracy than the spatial accuracy of 
GCMs — usually with good justification: spatial errors pose 
greater problems, especially for the short and medium range 
duration runs typically performed with the models. This ob- 
viously does not apply for long du ration runs, particularly i f 
quantitative predictions are sought ( |Thrastarson & Cho 2010). 

As already discussed, a computational mode arises in the 
leapfrog scheme, which is an example of a two-step scheme: 



xf/ 11 



= * n_1 + J"(* n ,x,f n ;e), 



(16) 



where xeR 2 ; recall that e is the Robert-Asselin time filter 
coefficient. Computational modes arise in all multistep meth- 




icr 10° 
Timestep 



FIG. 4. — Courant number as a function of time for two sets of runs with 
different values of r^, in each set. The four runs within each set have different 
e values, setting the strength of the Robert-Asselin time filter. For clarity each 
time series in a set has been offset vertically by 0.1 in the plot; and, the two 
sets, as groups, have been offset vertically by 0.5. The lower set of runs have 
To, = 3 Tp, while the upper set of runs have r t (, 
with e = 0.001 is indicated with a blue line, e - 
green line, and e = 0.01 a black line. 



= 0.1 Tp. For each set the run 
0.002 a red line, e = 0.006 a 



ods. Fortunately, in some multistep methods Af can be cho- 
sen to keep the amplitude of the modes from growing. How- 
ever, the leapfrog scheme is unstable for diffusion, for all Af . 
For this reason, the diffusion part of the equations is "time- 
lagged" by evaluating the diffusion terms at the time level 
(n- 1). This effectively time-marches the diffusion part by 
a first-order scheme. 

The Robert-Asselin filtered leapfrog scheme has been ana- 
lyzed by Durran ( 1991 1 for the simple oscillation equation, 



df" 



= iujip, 



(17) 



where lu is the frequency of oscillation. That analysis shows 
that, in the limit coAt <C 1, the relative phase-speed error of 
the physical mode is 



phys 



1 + 



l + 2e 
6(l-e) 



(wAf) 2 . 



(18) 



Therefore, the phase of the numerical solution leads the actual 
solution in time, and the error increases with larger e. No 
analysis exists to guide in choosing e. Hence, it is important 
to assess the sensitivity of the simulations to the filter. In this 
work, we have performed series of simulations in which e has 
been varied while keeping everything else fixed, for different 
values of relaxation time r t h. The simulation parameters are 
summarized in Table [3] 

Figure|4]shows the evolution of the Courant number [i* for 
simulations with different e, for two sets of runs with different 
values of r t h in each set. Note that for clarity each time series 
in a set has been offset vertically by 0. 1 in the plot; and, the 
two sets, as groups, have been offset vertically by 0.5. At 
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the T2 1 resolution of the runs shown, for t± = 3 t p and v = 
10 22 m 4 s" 1 , it is found that a value of at least e = 0.002 is 
needed to prevent the simulation from succumbing to time- 
splitting instability. 

With shorter relaxation time (t± = 0.1 t p ), a larger value of 
e is required for the simulation to proceed without blowing 
up; this is perhaps not surprising, in light of the preceding 
discussion. But, remarkably, even without explicit numeri- 
cal viscosity turned on (i.e., v set to 0), the simulation can 
proceed without crashing; and, this is so despite the fact that 
the physical field is completely swamped with noise! When 
T"th = 3 Tp, runs do not crash as long as e > 0.006. However, 
with r t h = 0.1 t p , the minimum e for not crashing is an order 
of magnitude greater. Evidently, an e value used in earlier 
studies of Earth's atmosphere should be adjusted when adapt- 
ing an Earth GCM for extrasolar planet study. In general, a 
shorter rth or lower viscosity requires stronger Robert- As selin 
filtering to prevent blow-up. 

Note that the Courant-Friedri chs-Lewy (CF L) criterion for 
stability of the leapfrog scheme (Durran[i999 1, 



< 



1 



can sometimes be exceeded in the middle of a run, even 
though the simulation is stable at t = (cf., run E3a in Fig- 
ure 0J. This is because the advective time-stepping limit de- 
pends on the maximum speed c max , which can increase during 
the evolution of a flow. Careful monitoring of the physical 
fields shows a zone of intense shear between the two vortices 
generating small-scale oscillations that rapidly amplify until 
the simulation becomes nonsense. The culprit is not lack of 
spatial resolution or a blocked turbulent cascade, in this case, 
because the calculation can be extended indefinitely by halv- 
ing the timestep. 

Figure [5] illustrates the sensitivity of the evolution to e for 
the range, e € [0.01,0.1]. Although the resolution in these 
calculations is only T21, they illustrate the point. Snapshots 
of the flow field are shown at three successive times for three 
simulations, differing only in the value of e. Similar flow pat- 
terns emerge in all the simulations: they all exhibit a cyclic 
behavior with vortices translating around the planet, undergo- 
ing large variations in strength and size as they do so, with 
corresponding changes in the temperature field. However, at 
a given instant the flow and temperature fields look different 
between the three runs. At t = 130r ; „ in all the runs there is a 
warm cyclone pair centered west of the substellar point. And 
in all the runs, the cyclones move westward and the flow and 
temperature fields undergo substantial changes before eventu- 
ally returning to a similar state, 15-20 planet rotations later. 
But at t = 145 T p , the run with the largest e has already re- 
turned to a state similar to that at t = 130r p , while the runs 
with smaller e take longer to complete their cycles. 

Figure [6] shows the behavior more clearly. The tempera- 
ture at a point on the model planet atmosphere (0° longitude, 
30° latitude) evolves in time for two simulations which have 
identical parameters, except for e. The two runs match nearly 
exactly until about 45 t p , when the two runs start to deviate. 
In the beginning only slightly. But, at about 70 t p the temper- 
ature oscillations in the run with the larger e lead in phase, 
compared to the run with smaller e. This behavior agrees 
qualitatively with Equation 18 Over long timescales the three 
simulations exhibit very similar behavior, even if amplitudes, 
phases, and periodicities of the flow and temperature fields 



are not exactly the same. As noted, simulations shown in Fig- 
ures [5] and [6] are at T21 resolution, but with higher resolution 

deviations appe ar even ear lier. 

In |Showman eTaLl ( [2009] , the MITgcm ( |Adcroft et al.|2004) 
is used. In that study, the model employs the third-order 
Adam s-Bashforth m ethod, which has some attractive prop- 
erties ( Durran|[l991[ l. However, the scheme does require an 
initialization phase in which ^ x and 'J 2 are computed from 
the initial condition by some other procedure, such as the 
fourth-order Runge-Kutta or a first- or second-order scheme 
with several short timesteps. It should be emphasized that — 
as it is a major point of this paper — the main concern is usu- 
ally adequate spatial resolution, especially in problems with 
inherent small-scale phenomena, not the time-integration. A 
second- or even first-order time-integration scheme can be 
perfectly adequate for many purposes. 

4. CONCLUSION 

A major aim of this paper has been to shed light on some 
crucial aspects for numerical modeling of atmospheric circu- 
lation on hot extrasolar planets. Here we have shown that, a 
spectral model offers advantages in accuracy and diagnostics, 
given that the higher-order field and wavenumbers are what's 
actually evolved. However, all numerical models, including 
spectral models, have limitations in how well they can rep- 
resent physical reality. Moreover, the models can easily be 
applied outside the realm of "safe parameters" and produce 
results that are nonsensical. The challenge is to properly test 
and identify the limits. When numerical artifacts appear, it is 
important to know how to deal with them and to know when 
a simulation should be discarded. 

In this paper we have shown that, for hot extrasolar planets 
simulations with stationary forcing, there is a strong sensitiv- 
ity to the strength of applied artificial viscosity. In addition, 
there is a relation between the thermal relaxation time r^, and 
the viscosity v. small r t h's lead to a large amount of unphysi- 
cal, grid-scale oscillations in the simulation, which forces the 
use of excessive amounts of artificial viscosity to quench the 
oscillations. Hence, using a fixed strength of artificial vis- 
cosity in a simulation with a large range of t± in the model 
domain (e.g., from about an hour to tens of days) — as done in 
many simulations in the literature — inevitably produces flow 
and temperature fields, which are either dominated by un- 
physical noise or over-damping. One may then wish to apply 
a spatially varying v, but clearly this is then motivated by a 
numerical basis rather than a physical one. 

The proper values to use for the relaxation time (or variables 
needed for realistic radiative transfer) are not known. Based 
on the findings in this work, calculations with extremely short 
T t h's warrant further scrutiny. Current GCMs may not be 
standing up too well to this stressful test. If, however, the 
short r t h are really physically relevant, then another form of 
heating/cooling parameterization or setup is needed. This is 
not a criticism of the Newtonian relaxation scheme, which in 
fact has been (and continues to be) very useful for understand- 
ing basic mechanisms. 

One solution could be to spatially vary the v, as already 
discussed; but, this would lead to further complexity. Even 
if direct radiative transfer is incorporated, one must ensure 
that the forcing is not too violent or strong (large amplitude 
and short timescale). Indeed, a comparison of our v values, 
scaled appropriately for the Earth, shows that we have had to 
use v values higher than that normally used in Earth studies 
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FIG. 5. — Temperature (color coded in K) with streamlines overlaid, for three simulations differing only in the value of e, shown at three moments in time. 
From left to right, the snapshots are taken at t = {130, 138, 145}t () . The top row is from a run with e = 0.01, the middle row with e = 0.06 and the bottom row e = 
0.10. All the fields are shown at the p S3 900 mb level. The substellar point is at 0° longitude and latitude. 
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FIG. 6. — Temperature at a fixed point, at 0° longitude and 30° latitude, as 
a function of time for the first 150 planet rotations in two simulations. The 
two curves show results of simulations that are identical apart only from the 
strength of the Robert- Asselin time filter. The red curve is from a run with e 
= 0.10, and the blue curve from a run with e = 0.01. 



( [Collins et aLl|2004) l. As discussed in ( |Chol[2008| ), if the ra- 
diative processes appear as practically instantaneous from the 
perspective of the flow, then an adiabatic approach is more ap- 
propriate. Certainly from a numerical accuracy standpoint, as 
motivated by the present work, adiabatic and "gently forced" 
calculations are useful as baselines. Else, gradually ramping 



up heating and/or initializing simulations close to a balanced 
state is necessary (Thrastarson & Cho|20101 l. 

GCMs of extrasolar planet atmospheres have great value in 
helping to guide and interpret observations. It is then impor- 
tant to critically examine the effects of the numerous param- 
eters that are specified. This is particularly crucial when ap- 
plying the models to a "new regime", where the physical con- 
ditions differ markedly from a traditional (e.g., Earth) one. In 
this paper we have shown examples of how a commonly-used 
forcing can steer GCMs to produce misleading results and 
how numerical expediencies, such as the Robert-Asselin fil- 
ter, can produce slewing frequ ency as well as the well- known 
damping and phase-errors ( |Durran|1991 Williams 2009]). In 
addition, we have discussed diagnostics procedures to better 
assess the quality of a simulation using the vorticity field and 
energy spectra. Reliance on spatial and temporal averages can 
effectively conceal telltale signs that a simulation is not trust- 
worthy. 

A simulation which is properly resolving the flow should 
approximately conserve energy for a long time. This should 
be so even if this property is not explicitly built into the dis- 
cretization algorithm as in the scheme of Arakawa ( 1966| [this 
scheme conserves the domain-integrated energy and enstro- 
phy (jC 2 ) m me nonlinear advection term]. For only then can 
we trust that a calculation is not artificially driven to an un- 
physical region in the solution space. 
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